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Abstract 

We present a biophysical approach for the coupling of neural network ac- 
tivity as resulting from proper dipole currents of cortical pyramidal neurons 
to the electric field in extracellular fluid. Starting from a reduced three- 
compartment model of a single pyramidal neuron, we derive an observation 
model for dendritic dipole currents in extracellular space and thereby for the 
dendritic field potential that contributes to the local field potential of a neu- 
ral population. This work aligns and satisfies the widespread dipole assump- 
tion that is motivated by the "open-field" configuration of the dendritic field 
potential around cortical pyramidal cells. Our reduced three-compartment 
scheme allows to derive networks of leaky integrate-and-fire models, which 
facilitates comparison with existing neural network and observation models. 
In particular, by means of numerical simulations we compare our approach 
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with an ad hoc model by Mazzoni et al. [Mazzoni, A., S. Panzeri, N. K. 
Logothetis, and N. Brunei (2008). Encoding of naturalistic stimuli by local 
field potential spectra in networks of excitatory and inhibitory neurons. PLoS 
Computational Biology 4 (12), el000239], and conclude that our biophysically 
motivated approach yields substantial improvement. 

Keywords: biophysics, neural networks, leaky integrate-and-fire neuron, 
current dipoles, extracellular medium, field potentials 



1. Introduction 



Sinc e Hans Berge r's 1924 discovery of the human electroencephalogram 
( EEG ) f lBergerl 119291 ) , neuroscientists achieve d much progress in clarifying 
its neural generators flCreutzfeldt et al.lll966al lbl iNunez and Srinivasanll2006 , 
Schomer and Lopes da Silval l201lh . These are the cortical pyramidal neu- 
rons, as sketched in Fig. [TJ that possess a long dendritic trunk separat- 
ing mainly excitatory synapses at the apical dendritic tree from mainly in- 
hibitory synapses at the soma and at th e perisomatic basal dendritic tree 
( iCreutzfeldt et al.lll966al . ISprustonl 120081 1 . In addition, they exhibit an ax- 
ial symmetry and are aligned in parallel to each other, perpendicular to the 
cortex' surface, thus forming a palisade of cell bodies and dendritic trunks. 
When both kinds of synapses are simultaneously active, inhibitory synapses 
generate current sources and excitatory synapses current sinks in extracellu- 
lar space, hence causing the pyramidal cell to behave as a microscopic dipole 
surrounded by its characteristic electrical field, the dendritic field potential 
(DFP). The densely packed pyramidal cells form then a dipole layer whose 
superimposed currents give rise to t he local field poten t ial (LFP) of neural 



masses and eventually to the EEG (ILinden et al .1120101. Linden et al.l 12011 



Nunez and Srinivasanll2006l . ISchomer and Lopes da SilvaJ 120111 ) . 



2 



apical dendritic tree 



perisomatic dendritic tree 




Figure 1: Sketch of a cortical pyramidal neuron with extracellular current dipole between 
spatially separated excitatory (open bullet) and inhibitory synapses (filled bullet). Neural 
in- and outputs are indicated by the jagged arrows. Dendritic current J D causes dendritic 
field potential (DFP). 



Despite of the progress from experimental neuroscience, theoretically un- 
derstanding the coupling of complex neural network dynamics to the electro- 
magnetic field in the extracellular space po ses challenging probl e ms; some of 

them have b een addressed to some e xtent by lBedard et al.l (120041 ). IBedard and Destexhe 
fl2009h . and IBedard and Destexh^ fl2012h . 

In computer simulation studies, neural mass potentials, such as LFP 
and EEG are most realistically simulated by means of multicompartmen- 



ts 

tal models (ILinden et al.l 12010. ILinden et al.l 120111 . iProtopapas et al.l 11998 



Sargsvan et aDl200lf ). ILinden et all fl2010f ) calculated the current dipole mo 



mentum of the DFP for single pyramidal and stellate cells, based on several 
hundreds compartments of the dendritic trees. Their results were in com- 
pliance with the standard dipole approximation of the electrostatic multi- 
pole expansion in the far-field (more than 1 mm remote from the dendritic 
trunk), but they found rather poor agreement with that approximation in 
the vicinity of the cell body. For comparison they also computed a "two- 
mo nopole" model of one synaptic current and its counterpart, the somatic 
return current, estimated from the current dipole momentum of the whole 
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dendritic tree. This "two-monopole" model, which corresponds to an electri- 
cally equivalent single dipole model, obtained from the decomposition of the 
dendrite into two compartments, better approximates the true current dipole 
momentum in the vicinity of the pyramidal neuron. By superim posing; the 
DFPs of pyramidal cells to the ensemble LFP. lLinden et al.l (120 llh found that 
LFP properties cannot be attributed to the far-field dipole approximation. 

However, realistic multicompartmental models are computationally too 
expensive for large-scale neural network simulations. Therefore, various 
techniques have been proposed and employed to overcome computational 
complexity. These include networks of point models ( i.e. devoid from any 
spatial representation), b ased on conductance models (IHodgkin and Huxley 
1952 . Mazzoni et al. 20081), population dens i ty mo dels ( Omurtag et al. 2000 ). 
or firing rate models (jWilson and Cowan! Il972l ). which can be seen as a 
sub class of populat ion density models, with uniform density distribution 
(IChizhov et al.ll2007l ). In these kinds of models, mass potentials such as LFP 
or EEG are conventionally described as averaged membrane potential. A 



different class of models are neural mass models ( 


David and Friston 


2003. 


Jansen and Rit 


1995. 


Wendline et al. 


2000, 


Rodrieues et al. 


20 id), where 



mass potentials are esti mated either through sum s (or actually differences) of 
postsynaptic potentia ls ( David and Friston 20031 ) or of postsynaptic currents 



(IMazzoni et al.l 12008 



In particular, the model of lMazzoni et al.l ( 12008) which is based onlBrunel and Wang 



( 120031 ). recently led to a series of follow-up studies ( IMazzoni et al.ll2010l . l20 111 ) 
addressing the correlations between numerically simulated and experimen- 
tally measured LFP/EEG with spike rates by means of statistical modeling 
and information theoretic measures. In all of the above point models and 
their extension to population models, it is assumed that the extracellular 
space is iso-potential and the majority of studies thereby neglect the effect 
of extracellular resistance. That is, the extracellular space constitutes a dif- 
ferent and isolated domain with no effect on neur onal dynamics. 

In this article we extend the ad hoc model of lMazzoni et al.l ( 120081 ) to- 
wards a biophysically better justified approach, taking the dipole character 
of extracellular currents and fields into account. Basically, our model corre- 
spond s to the "two-monopole", or, equivalent dipole model of iLinden et al 



( I2010l ) which gave a good fit of the DFP close to the cell body of a cor- 
tical pyramida l neur on. However, we aim to keep the simplicity of the 



Mazzoni et al.l ( 120081 ) model in terms of computational complexity, by en- 



dowing the extracellular space with resistance and by keeping point-like neu- 
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ronal circuits. That is, in our case we do not quite consider point neurons, 
nor spatially extended models with detailed compartmental morphology, yet 
an intermediate level of description is achieved. To this end we propose a 
reduced three- compartmental model of a si ngle pyramidal neuron (IDestexhe 
200ll . lbeim Grabenll2008l . IWang et al.ll2004l ). and derive an observation model 
for the dendritic dipole currents in the extracellular space and thereby for the 
DFP that contributes to the LFP of a neural population. Interestingly, our 
reduced three-compartmental model enabl es us to derive a lea ky integrate- 
and-fire mechanism (as for a point model (IMazzoni et al.ll2008l )). with addi- 
tional observation equations for the DFP, which all together allows to study 
the relationship between spike rates and LFP. Our derivations also nicely map 
re alistic electrotqnic pa rameters to phenomenological parameters considered 



in 



Mazzoni et all ( 120081 ) . 



2. Material and Methods 



Mazzoni et al.l (120081 ) consider three populations of neurons, namely exci- 
tatory cortical pyramidal cells (population 1), inhibitory cortical interneurons 
(population 2) and excitatory thalamic relay neurons (population 3), pass- 
ing sensory input to the cortex that is simulated by a random (Erdos-Renyi) 
graph of K = 4000 pyramidal and L = 1000 interneurons with connection 
probability P = 0.2. 



2.1. Theory 

We describe the ith cortical pyramidal neuron [Fig. [T] from population 
1 via the electronic equivalent (reduced) three-com partment model Fig. [2] 
(|Pestexhell20"0ll . Ibeim Graben1l2008l . IWang et al.ll2004 ) , which is parsimonious 
to derive our observation model: one compartment for the ap ical dendritic 



tree, another one for soma and perisomatic basal dendritic tree (ILinden et al 
20101 ) . and the third — actually a leaky integrate-and-fire (LIF) unit — for 



the axon hillock where membrane potential is converted into spike trains by 
means of an integrate-and-fire mechanism. 
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Figure 2: Proposed electronic equivalent circuit for a pyramidal neuron (reduced three 
compartmcntal model). Note that the apical and basal dendrites are not true compart- 
ments since capacitors are not explicitly represented, rather, these are implicitly taken 
into account via EPSP and IPSP static functions, thus keeping computational complexity 
low. 



Excitatory synapses are represented by the left-most branch, where exci- 
tatory postsynaptic potentials (EPSP) at a synapse between a neuron j from 
population 1 or 3 and neuron i act as electromotoric forces Efy These poten- 
tials drive excitatory postsynaptic currents (EPSC) If^ essentially consisting 
of sodium ions, through the cell plasma with resistance Rfj from the synapse 
towards the axon hillock. 

The middle branch describes the inhibitory synapses between a neuron 
k from population 2 and neuron %. Here, inhibitory postsynaptic potentials 
(IPSP) E] k provide a shortcut between the excitatory branch and the trigger 
zone, where inhibitory postsynaptic currents (IPSC) l\ k (essentially chloride 
ions) close the loop between the apical and perisomatic dendritic trees. The 
resistivity of the current paths along the cell plasma is given by R\ k . 

The cell membrane at the axon hillock itself is represented by the branch 
at the right hand side. Here, a capacitor Cj reflects the temporary storage 
capacity of the membrane. The serial circuit consisting of a battery E M 
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and a resistor R M denotes the Nernst re sting potential and the l eakage con 



ductance of the membrane, respectively (IJqhnston and Wu 



19971) . Finally, a 



spike generator flHodgkin and Huxleyl 119521 . iMazzoni et al.l 120081 ) (indicated 
by a "black box") is regarded of having infinite input impedance. Both, 
EPSP and IPSP result from the interaction of postsynaptic receptor kinetics 
with dendritic low-pass filtering in comp artments one and two, respectively 



(IDestexhe et al.lll998l . iLinden et al.ll2010l ). Hence the required capacitances, 
omitted in Fig. El are already taken into account by E®, E\ k . Therefore, we 
refer to our model as to a "reduced compartment model" here. 

The three compartments are coupled through longitudinal resistors, Rf , Rf, Rf, Rf 
where R^,Rf deno te the resistivity of t he cell plasma and Rf,Rf that of 
extracellular space (IHolt and Kochlll999l ). 

Finally, the membrane voltage at the axon hillock Ui (the dynamical state 
variable) and the DFP Vi, which measures the drop in electrical potential 
along the extracellular resistor Rf are indicated. For the aim of calculation, 
the mesh currents If (the dendritic current), If (the basal current) and lj F 
(the integrate-and-fire current) are indicated. 

The circuit in Fig. [2] obeys the following equations: 



IF 



17 



E 



E ik 



/Ah 



k=l 



c, 



dUi Ui-E 



dt 



R 



Rpfi + W + + (Rf + + u i ,i<j< P 



13 V 
>I rl 



R\ k l\ k + (Rf + Rf)lf +Ui ,1 < k < q 
R?I? . 



(1) 

(2) 

(3) 
(4) 

(5) 
(6) 
(7) 



Here, p is the number of excitatory and q is the number of inhibitory synapses 
connected to neuron i. 

The circuit described by Eqs. (p]-[7]) shows that the neuron i is likely to 
fire when the excitatory synapses are activated. Then, the integrate-and-fire 
current lj F equals the dendritic current if. If, by contrast, also the inhibitory 
synapses are active, the dendritic current if is shunted between the apical 
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and perisomatic basal dendritic trees and only a portion could evoke spikes 
at the trigger zone [Eq. fll}]. On the other hand, the large dendritic current 
If flowing through the extracellular space of resistance Rf, gives rise to a 
large DFP V t . 

In order to simplify the following derivations, we gauge the resting po- 
tential [Eq. (3])] to E M = 0, yielding 



From we obtain the individual EPSC's as 
And accordingly, the individual IPSC's from (jfJJ) 



Inserting (J9j) into ([I]) yields the excitatory dendritic current 

= E - + W + (i? + R?)i? + m , (11) 

3=1 V 

where we have introduced the excitatory dendritic conductivity 

sf = tw- (12) 

Likewise we obtain the inhibitory dendritic currents from ([2]) and (fTOj) as 

J * = E TT^ - ^t ( ^ B + ^ C)J ^ F + ^ ' (13) 

fc=l n ik 

with the inhibitory dendritic conductivity 



k=l 



s 



With these results, we obtain an interface equation for an observation 
model as follows. Rearranging ( II ip yields 



+ sfW + ^ D )] = E -Le% - gf[(Rf + R?)if + m (is) 
i=i ^' 

Next, we eliminate lj F through (JH]): 

If [1 + (R? + fl?)] = £ ^ " sf h (iff + fl, c ) ^ + 14 (l + £±£) 

j = l l 3 



Division by 1 + gf (Rf + R?) gives the desired expression for the extra- 
cellular dendritic dipole current: 



A D = E«,,^-A^-7^, (16) 

3=1 

with the following electrotonic parameters 

aid = W+WW^Rf)] (17) 

n Ci9?(Rf + ffi) f-,o\ 

Pi l + gf(Rt + R?) 1 J 

7i ^[i + ^ + flO)]' ^ 

In order to derive the evolution equation we consider the integrate-and- 
fire current lj F that is given through (J3j) . The individual EPSCs and IPSCs 
have already been obtained in (jHj) and (TTDjl . respectively. Inserting ( JT51) into 
(ED yields 



Next we insert our interface equation (fTB]) and also 



1 dt + i? M 



[l-gXRf + Rfy-glUt = j^iSg-fl 

3=1 



dC^ 
dt 



9 1 
Ui ~ E 



pi 

k=l U ik 
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and obtain after some rearrangements 



{Q [1 - gj(R? + Rf)] + Pi} 



dt 



+ 



1 - g](Rf + Rf + R M ) + R M j 



R M 
q 



1 



3=1 



pi 

k=l 



and after multiplication with 



R 



M 



1 1 - g}(Rf + Rf + R M ) + R M 7i 

the dynamical law for the membrane potential at axon hillock: 

p <? 



1 dt 



U i = ^2wf j Ef j -^2w] k E} k , 



k=l 



where we have introduced the following parameters: 

• time constants 

r l = n {C l [l-g\{Rf + Rf)\+(5 l } 

• excitatory synaptic weights 



inhibitory synaptic weights 



w 



ik 



R \k 



(20) 



(21) 



(22) 



(23) 



Using the result (TSUI) , we can also eliminate the temporal derivative in 
the interface equation ffl6|) through 



dt T: 



-17, 

Li=i 



k=l 



(24) 



which yields 



» E 



E A 7~iE X ^ ^ I 7~>I ft' 

' k=\ Ti 



T, 



li)Ui. 
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And eventually, by virtue of Eq. (0) after multiplication with Rf the dendritic 
field potential (DFP) 



v t = j2 E% + E ^ + fc* 7 * ' ( 25 ) 

j=i fc=i 



with parameters 



tig = flj>«J(I-£) (26) 



<4 = «, D «4^ (27) 
£ = f- - 7^ • (28) 



The change in sign of the inhibitory contribution from Eq. (|20|) to Eq. (|25|) has 
an obvious physical interpretation: In ( 120]) . the change of membrane potential 
Ui and therefore the spike rate is enhanced by EPSPs but diminished by 
IPSPs. On the other hand, the dendritic shunting current if in (j25p is large 
for both, large EPSPs and large IPSPs. 

From Eq. ( 120]) we eventually obtain the neural network's dynamics by 
taking into account that postsynaptic potentials are obtained from presy- 
naptic spike trains through temporal convolution with postsynaptic impulse 
response functions, i.e. 



Ef{t)= f sf\t-t')R 3 {t')dt' 

J — oo 



(29) 



Ell 

where (£) are excitatory and inhibitory synaptic impulse response func- 
tions, respectively, and Rj is the spike train 

R j {t) = '%26(t-t v -TL) (30) 

u 

coming from presynaptic neuron j, when spikes were emitted at times t u . 
The additional time c onstant tl is attributed to synaptic transmission delay 
( IMazzoni et al.l 120081 ). These events are obtained by integrating (120]) with 
initial condition 

Ui{t v ) = E. (31) 
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Where E is some steady-state potential ( jMazzoni et al.l 120081 ) . If at time 
t = t,. the membrane reaches a threshold 



(32) 



Ui(t) > 9i(t) 

(with possibly a time dependent activation threshold 9i(t)) from below 
then an output spike 8(t — t v ) is generated, which is then followed by a 
potential resetting as follows 



dUjjt) 
dt 



> 



+1, 



<- E 



(33) 



Additionally, the integration of the dynamical law is restarted at time t = 
t u+ i + r rp after interrupting the dynamics for a refractory period r rp . 

Inserting ff29l) into f|20l) entails the evolution equation of the neural net 
work 

dU, 
dt 



(34) 



fc=i 



where the signs had been absorbed by the synaptic weights, such that wfj > 



for excitatory synapses and w\ k < for inhibitory synapses, respectively 

Ell 

Following iMazzoni et al.l (120081 ) an individual postsynaptic current 1^ 
at a synapse between neurons i and j obeys 



T. 



E| 



id/, 



Ell 



E|l 



E|l 



dx 



Ell 



T, 



E|I ^jj , E|I 
dt + ^' 



.r 



F E|1 



(35) 
(36) 



where rj' 1 are decay time constants and r^' 1 are rise time constants of EPSC 
and IPSC, respectively. Auxiliary variables are denoted by xfj 1 , while F®' 1 
prescribes presynaptic forcing 



Ff = nJ^it) 



(37) 



Ell 

v denotes synaptic gain with v = 1 mV 



with spike train fl3"01) . Here, Jij 
as voltage unit. 

Note that (|37|) is essentially a weighted sum of delta functions, such that 
a single spike can be assumed as particular forcing 



F S{t) 



(38) 
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with some constant F . 

Derivating (j35p and eliminating xfj 1 transforms Eqs. (1331 - 1331) into a 
linear second-order differential equation with constant coefficients 

rf'rf^ + (rf + rf^ + if = F f . (39) 

Equation (1391 with the particular forcing (1381) is solved by a Green's 
function sf l (t) such that the general solution of (1391) is obtained as the 
temporal convolution 

/f (*) = f sf\t-t')Ff{t)dt' . (40) 

»/ — oo 

For t 7^ 0, (1391) assumes its homogeneous form and is easily solved by 
means of the associated characteristic polynomial 

^ E|I A 2 + (rf + r r E l I )A + l = (41) 

with roots Ai = — 1/rJ' 1 and A2 = — 1/Tr, entailing the Green's functions 

sf\t) = (a^ 1 ^ 1 - B^ 1 ^) e(t) (42) 

with the Heaviside step function 0(t). 

The constants A E \\B E \ l > are obtained from the initial conditions 

Ell 

s i (t) — 0, reflecting causality, and a suitable normalization 

00 

sf\t)dt = 1 . 







The initial condition yields A E ^ = B E ^ 1 = S' E ' 1 , while the remaining 
constant 

oE|I = 1 

u Ell Ell ' 

d ue to normalization . Ther efore, the normalized Green's functions are those 
of iBrunel and Wangl (120031 ) 

fW = ^(e^-e^)e(t). (43) 
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Now, we are a ble to compare our DFP Vi [Eq. (1251) ] with the estimate 
of iMazzoni et al.l (120081 ) which is given (in our notation) as the sums of the 
moduli of excitatory and inhibitory synaptic currents, i.e. 



Vi 



3 LB _ | ; 

j 



k 



(44) 



where "MPLB" refers to the authors IMazzoni et al.l (120081 ). 

From (l2l)|) and (l4"4l) . respectively, we compute two models of the local 
field potential (LFP). First, by summing DFP acros s all pyramidal neurons 
(Ibeim Graben and Kurthsll2008l. IMazzoni et al.ll2008r) . and, second by taking 
the DFP average ( Nunez and Srinivasanll2o"o6 ). which yields 



In 


_ ^MPLB 

i 


(45) 


Li 


1 taMPLB 

K * 

i 


(46) 


L 3 




(47) 


L 4 




(48) 



where K is number of pyramidal neurons. 

2.2. Parameter estimation 

Next, we relate the elect rotonic parame t ers of our model to the phe- 
nomenological parameters of IMazzoni et al.l (120081 ) . To this end, we first 
report their synaptic efficacies in Tab. [TJ 

Synaptic efficacies / mV on interneurons on pyramidal neurons 



GABA 2.7 
recurrent cortical AMPA 0.7 
external thalamic AMPA 0.95 



1.7 
0.42 
0.55 



Table 1: Parameters laid as in IMazzoni et al.l (|2008l ) 



From these, we compute the synaptic weights through 



0.42 if j "cortical" 
0.55 if j "thalamic" 



(49) 
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and 



w 



ik 



Ji/v = 1.7 



Next, we determine the factors r; by virtue of Eq. ([23]) through 



W; 



ik 



1.7 



S'gaba 1 nS 



1.7GQ 



using the inhibitory synaptic conductivity (?gaba = 1 nS, Correspondingly, 
Eq. (122]) allows us to express ajj in terms of the excitatory synaptic weights 
through 



w 



a 



•J 



0.25 nS if j "cortical" 
0.32 nS if j "thalamic" 



From aij we can determine the total excitatory synaptic conductivities 
gf according to Eq. (TTTI) through 



a 



R?)} 



i - (Rf + R?) E 



'j 



and hence 



1 



Oy[l+^(i^ + i2?)] 

Inserting next ffl8|) into (T2T]) yields 

1 + gf (itf + i??) + (i2f + R?){gf - g][l + gf{Rf + Rf))\ 



1 + gf{Rf + R\ 



(50) 
(51) 

(52) 



Equation (1521) could co nstraint the choice o f the membrane capacitance 
Ci by choosing = 20 ms (iMazzoni et al.ll2008l ). 

In order to also determine the DFP parameters (126]) - (128]) . we finally 
compute the ratios 



Pi 



g?(R? + R?) 



n{i + gf(Rf + R?) + (Rf + R?){gf - + g?{Rt + R?)}}} 
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The remaining electrotonic parameters Rf 1 , Rf-,Rf,Rf, and Rf are es- 
timated from cell geometries as follows. The resistance R of a volume con- 
ductor is proportional to its length i and reciprocally proportional to its 
cross-section A, i.e. 

R = P j (53) 

where p is the (specific) resistivity of the medium. Table [2] shows the resistiv- 
ities of the three kinds of interest which then allows to evaluate the various 
volume conductor resistances according to Eq. (1531) . 



medium 



p/Qcm 



cell membrane (at axon hillock) 5 • 10 7 
cell plasma (cytoplasm) 200 
extracellular space 333 

Tabl e 2: Resistivities of cell membrane, c ell plasma and extra c ellular space . Para meters 
from lKole and Stuartl (|2012f) . iRalll (|l977f ). iMainen et all (fl995h . iGold et all (j2007l) . Note 
that the resistivity of the cell membrane has to be related to the constant membrane 
thickness (sa 10 nm). 



We consider a total dendritic length of 21 = 20 /im and a dendritic radius 
of a = 7/im, that are generally subjected to variation. Equally, parameters 
that were allowed to vary are the length and radius of the axon hillock, 
yet herein we consider a length of 21 = 20 ixm and radi us of a = 0.5 /im 
f Destexhe 200 ll . Kole and Stuart 2012 . Mainen et al. 1995 ). To evaluate the 
intracellular (Ra, Rb) and extracellular (Rd, Rc) resistances, respectively, 
according to Eq. (153]) . we consider a simple implementation where the length 
i is half of the dendritic length (i.e. basal and apical length are symmetrical, 
but this can be broken). However, the cross sectional area for the cytoplasm 
is simply A = na 2 . Finally, the area of the axon hillock is simply the surface 
area of a cylinder. 

In order to also determine the cross-section of extracellular space between 
dendritic trunks we make the following approximations. We assume that 
dendritic trunks are parallel aligned cylinders of radius a and length I that 
are hexagonally dense packed. Then the centers of three adjacent trunks 
form an equilateral triangle with side length 2a and hence area 2 v / 3a 2 . The 
enclosed space is then given by the difference of the triangle area and the 
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area of three sixth circle sectors, therefore 



A 



space 



-ira 
6 



Hence, the cross- section of extracellular space surrounding one trunk is 

(54) 



A = QA spacc = (l2y/3 - 3tt 



2. 3. Simulations 



Subsequently, we imp lement an identical network to the one considered 
by iMazzoni et al.l (120081) with Brian S imulator, that is a Python based en- 
vironment ( IGoodman and Brettell2009l ). However, the derivations from the 
previous section enables the possibility of setting a dipole observable that 
measures the local dendritic field potential (DFP) on each pyramidal neurons, 
given by Eq. ( 1251) . This allows then to define a mesoscopic LFP observable, 
which can be equated either as averaged DFP or simply given as the sum of 
DFP, given by Eqs. (1451 - 148]) . Primarily, we compare our LFP measure L4, 
proposed as the average of DFP, with the IMazzoni et al.l LFP L\ which is 
defined as the sum of absolute values of GABA and AMPA currents Eq. (1441) . 
Additionally, we al so compare all p ossible measures, na mely, mean mem brane 



potential Uj, IMazzoni et al.l LFP L±, average of IMazzoni et al.l DFP 

L2, sum of DFP L3, and the average of DFP L4. 

For completeness, we briefly summarize the description of the network 
(we refer the reader to IMazzoni et al.l (120081 ) for details) . The network mod- 
els a cortical tissue with leaky integrate-and-fire neurons, composed of 1000 
inhibitory interneurons and 4000 pyramidal neurons, which are described by 
the evolution equation ( 134|) . The threshold crossings given by Eq. (|32|) is 
considered static with 6*.; = 18 mV and the reset potential E = 11 mV. The 
refractory period for excitatory neurons is r rp = 2 ms while for inhibitory 
neurons it is r rp = 1 ms. The network connectivity is random and sparse 
with a 0.2 probability of directed connection between any pair of neurons. 
The evolution of synaptic currents, fast GABA (inhibitory) and AMPA (ex- 
citatory) are described via the second order evolution equations (1351 - 136|) . 
which are activated by incoming presynaptic spikes represented by Eq. fl30|) . 
The latency of the postsynaptic currents is set to tl = 1 ms and the rise and 
decay times are given by Tab. [31 
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Synaptic times r r / ms / ms 

GABA 025 5 

AMPA on interneurons 0.2 1 

AMPA on pyramidal neurons 0.4 2 



Table 3: Synaptic rise (r r ) and decay times (r<j). Parameters laid as in iMazzoni et al 



((20081) 



Moreover, synaptic efficacies, , f° r simulation were presented in Tab. [I] 
Note that relation (1491) then allows to determine the synaptic weights. Addi- 
tionally, all neurons receive external thalamic excitatory inputs, that is, via 
AMPA-type synapses, which are activated by random Poisson spike trains, 
with a time varying rate that is identical for all neurons. Specifically, the 
thalamic inputs are the only source of noise, which attempts to account for 
both cortical heterogeneity and spontaneous activity. This is achieved by 
modeling a two level noise, where the first level is an Ornstein-Uhlenbeck 
process superimposed with a constant signal and the second level is a time 
varying inhomogeneous Poisson process. Thus, we have the following time 
varying rate, A(t), that feeds into inhomogeneous Poisson process: 

Tn ^d^ = - n ( t ) + a n^v(t) (55) 

X(t) = [co + n(*)]+ (56) 

where r)(t) represents Gaussian white noise, cq represents a constant signal 
(but equally could be periodic or other), and the operation [■] is the threshold- 
linear function, [x] + = x if x > 0, [x] + = otherwise, which circumvents 
negative rates. The constant signal cq can range between 1.2 to 2.6 spikes/ms. 
The parameters of the Ornstein-Uhlenbeck process are r n = 16 ms and the 
standard deviation a n = 0.4 spikes/ms. 

For complete exposition, we note that from an implementation viewpoint 
(within the Brian simulator), a copy of the postsynaptic impulse response 
function ( 129]) has to be evaluated to calculate the dendritic field poten- 
tial (DFP) (I25p with weights wfj- This implies evaluating the second or- 
der process ( 135]) - ( 136]) with a different forcing term. Specifically, starting 
from lfj(t) = wfjEfj 1 ^) = sf I (t) * F^ 1 and pre- multiplying both sides 

Ell 

with and subsequently re-arranging we obtain the desired forcing term 



F^ 1 = wfjF^/wfJ 1 . Note further that by expanding the term F^ 1 with 
equation ( 157|) and using relation (|4"9"]) we finally obtain Pf^ 1 = w^ l TivRj(t). 



18 



3. Results 



Following iMazzoni et al.l (120081 ). the network simulations are run for two 
seconds with three different noise levels, specifically, receiving a constant 
signal with three different rates 1.2, 1.6 and 2.4 spikes/ms as depicted in 
Fig. [3J Note that these input rates do not mean that a single neuron fires 
at these high rates. Rather, it can be obtained from multiple neurons that 
jointly fire with slower, yet desynchronized, rates converging at the same 
postsynaptic cell. The Poisson process ensures that this is well represented. 
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Signal rate = 1.2 spikes/ms 



Signal rate = 1.6 spikes/ms 



Signal rate = 2.4 spikes/ms 
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Figure 3: Dynamics of the network and LFP comparisons: The three columns represent 
different runs of the network for three different rates, 1.2, 1.6 and 2.4 spikes/ms. In each 
column, all panels show the same 250 ms (extracted from 2 seconds simulations). The top 
panels (A-C) represent thalamic inputs with the different rates. The second top panels 
(D-F) corresponds to a raster plot of the activity of 200 pyramidal neurons. Panels (G-I) 
depict average instantaneous firing rate (computed on a 1 ms bin) of interneurons (blue) 
and panels (J-L) corres pond to average instantaneous firing rate of pyramidal neurons. 



Panels (M-O) show the iMazzoni et all LFP L\ from Eq. ([4"5]l. Finally, panels (P-R) 
depict our proposed LFP measure L4, which is the average of dendritic field potential 
(DFP) [Eq. 
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The focus is to compare our proposed measure L4, defined as me an of 



the dendritic field potential (DFP) [Eq. (HB)) ]. with the Mazzoni et all LFP 



L\ from Eq. (145 p . In Fig. [3] one sees two main striking differences between 
the two measures, namely in frequency and in amplitude. Specifically, L\ 
responds instantaneously to the spiking network activity by means of high 
frequency oscillations. Moreover, Li also exhibits a large amplitude. In 
contrast, our mean dendritic field potential L 4 measures comparably to ex- 
perimental LFP, that is, in the order of millivolts, and although it responds 
to population activity, it has a relatively smoother response. Actually one 
can realize that our LFP estimate represents low-pass filtered thalamic input. 

The phy siological relevance of this is not yet clear in our work. However, 
recent work (IPoulet et al.ll2012l ) shows that desynchronized cortical state dur- 
ing active behavior is driven by a centrally generated increase in thalamic 
action potential firing (i.e. thalamic firing controls cortical states). Thus, it 
seems that cortical synchronous activity is suppressed when thalamic input 
increases, thereby suggesting that cortical desynchronized states to be re- 
lated to sensory processing. This work further quantifies these observations 
by applying Fast Fourier Transform (FFT) to cortical EEG and subsequently 
comparing with thalamic firing rate by means of Pearson correlation coeffi- 
cient. Unfortunately they do not quantify the amount of thalamic oscillations 
contained within the cortical EEG. 

Yet, to keep a c omparable comp arison between measures, we also compute 
the average of the iMazzoni et al.l DFP L 2 [Eq. f!48p ] and additionally the 
mean membrane potential (the standard considered in the neuroscientific 
literature). These are shown in Fig. HI 



21 



Signal rate = 1.2 spikes/ms 



Signal rate = 1.6 spikes/ms 



Signal rate = 2.4 spikes/ms 
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Figure 4: Comparison of different LFP measures when the network receives constant signal 
with three different rates (1.2, 1.6 and 2.4 spikes/ms). Again, only 250 ms is represented 
(extracted from 2 sec simulation). The hrst plot (A-C) corresponding to the the different 
rates shows the most widespread LFP measure used in the literat ure, namely av erage 



membrane potential ^i' The second panel (D-F) shows the IMazzoni et al.l LFP 

L\ from Eq. (|45|). The third panel (G-I) displays the average of the IMazzoni et aL DFP 
L 2 [Eq. (|46|) ]. Similarly, the fourth panel (J-L) shows the total, L 3 , [Eq. (j47| ] and the 
last panel (M-O) depicts the averaged, L4, [Eq. (|48|) ] LFP measure. Note the different 
amplitude scales between measures. 
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Clearly, in terms of time profile, the summed and averaged observables 
are similar within the same class of LFP measures. However, in all cases 
the iMazzoni et al.l LFP L\ exhibits a significantly larger order of magnitude, 
which diverges substantially fr om experimental LFP amplitudes, typi cally 
varying between 0.5 to 2 mV (jNiedermeyerl 12005k lLakatos et al.ll2005l ). In 
contrast, although the mean DFP is not contained within the interval from 
0.5 to 2 mV it arguably performs better. However, we do concede further 
work is required. Some gains in improving the different LFP measures can 
be achieved by applying for example a weighted average, which would mimic 
the dis tance of an electrode to a p articular neuron by means of a lead field 
kernel ( Nunez and Srinivasanl 2006 ). For example a convolution of either L\ 
or L 2 with a Gaussian kernel (representing the distance to a neuron), would 
yield a measure that captures better the local field potential or better the 
dendritic field potential of the nearest neurons. However, further work will 
be required to properly quantify the gain when space is taken into account. 

In Fig. [5] we finally contrast the power spectra of the different LFP mea- 
sures. 



23 




Figure 5: Comparison of power spectra of the various LFP measures when the network 
receives constant signal with three different rates (1.2, 1.6 and 2.4 spikes/ms): The first 
plot (A-C) corresponding to the different rates shows the power spectrum of the average 
membrane potential If J^i^i- The second plot (D-F) and third plots (G-I) show powe r 



spectra of the total and average of L\ and L2 corresponding to IMazzoni et al.l (|2008l ). 



respectively. The fourth plot (J-L) and fifth plots (M-O) display power spectra of the L3 
and L4 measures from our model, respectively. Note we show the full spectrum up to 
5 kHz only for convenience due to the fine sample rate. 



One interesting feature is that the power spectrum of the IMazzoni et al 



LFP measures decays much more slowly that the average membrane potential 
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for higher frequencies. This observation is true for both, Li and L 2 . In 
contrast, our LFP measures L 3 and L 4 fare better, and in particular, L 4 
decays at an approximately similar rate as the average membrane potential. 



4. Discussion 



In this article we derived a model for cortical dipole fields, such as den- 
dritic and local field potential (DFP/LFP) from biophysical principles. To 
that aim we decomposed a cortical pyramidal cell, the putative generator 
of those potentials, into three compartments: the apical dendritic tree as 
the place of mainly excitatory (AMPA) synapses, the soma and the peri- 
somatic dendritic tree as the place of mainly inhibitory (GABA) synapses 
and the axon hillock as the place of wave-to-spike conversion by means of an 
integrate-and-fire mechanism. From Kirchhoff 's laws governing an electronic 
equivalent circuit of our model, we were then able to derive the evolution 
equation for neural network activity Eq. (|34|) and, in addition, an observa- 
tion equation Eq. ( 125]) for the dendritic dipole potential contributing to the 
LFP of a cortical population. 

In order to c ompare our approach with anot her model discussed in the 
recent literature ( Mazzoni et al]2008 . 2010 , 2011) w e aligned the parameters 
of our model with the model of iMazzoni et al.l ( 120081 ) who approximated DFP 
as the sum of moduli of excitatory and inhibitory synaptic currents Eq. (144]) . 
From b oth approaches , we computed four different LFP esti mates: L-\ , th e 



sum of IMazzoni et al.l DFP, L 2 , the population average of IMazzoni et al. 
DFP, L 3 the sum of our dipole DFP, and L 4 the population average of our 
dipole DFP [Eqs. (g5]-|48D]. 

Our resul ts indicate two main effects between our dipole LFP measures 



and th ose of IMazzoni et al.l Firstly, the measures based on IMazzoni et al 



( 120081 ) systematically overestimate LFP amplitude by almost one order of 
magnitude. One reason for that could be attributed to the direct conversion 
of synaptic current into voltage without taking extracellular conductivity into 
account, as properly done in our approach. Yet, another, even more crucial 
reason is disclosed by our equivalent circuit Fig. |2j In our approach there is 
just one extracellular current I D flowing from the p erisomatic to the apical 
dendritic tree. In the model of IMazzoni et al.l ( 120081 ). however, two synaptic 
currents that might be of the same order o f magnitude ar e supe rimposed 
to the DFP. Secondly, the measures based on IMazzoni et al.l ( 120081 ) also sys- 
tematically overestimate LFP frequencies. This could probably be attributed 
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partly to spurious higher harmonics introduced by computin g absolute val- 
ues. Moreover, taking the power spectrum shows that the iMazzoni et al. 



(120081 ) measure decays much more slowly than the average membrane poten- 
tial, which is at variance with experimental data. 

However, at the current stage, both models, that of IMazzoni et al.l (|2008[ ) 
and our own, agree with resp e ct to the polarity of DFP and LFP. The mea- 
sures based on IMazzoni et al.l ( 120081 ) have positive polarity simply due to the 
moduli. On the other hand, also the direction of current dipoles in our model 
is constrained by the construction of the equivalent circuit Fig. [2] where cur- 
rent sources are situated at the perisomatic and current sinks are situated 
at apical dendrite tree. Taking this polarity as positive also entails positive 
DFP and LFP that could only change in strength. However, it is well known 
from brain anatomy that pyramidal cells appear in at least two layers, III 
and VI, of neocortex. This is reflected in experiments when an electrode tra- 
verses different layers by LFP polarity reversals, and, of course, by the fact 
that LFP and EEG oscillate between positive and negative polarity. Adapt- 
ing our model to this situation could be straightforwardly accomplished in 
the framework of neural fi eld theory by fully representing space and simulat- 



ing l ayered neural fields ( lAmaril 119771 iJirsa and Hakenl Il996l . Ibeim Graben 



20081). By contra s t such a generalization is impossible at all with the model 



of IMazzoni et al.l ( 120081 ) due to the presence of absolute values. 

On theses grou nds we have go od indication that our measure is an im- 
provement to the IMazzoni et al.l LFP measures, and, quit e importantly 



i t is b iophysically better motivated than the ad hoc model of IMazzoni et al 
(120081 ). However, much considerable effort is still required to underpin all the 
relevant LFP mechanisms and to better represent experimental LFP/EEG 
dynamics. 

Finally, our work provides a new framework where dendritic field poten- 
tials and the relationship between firing rates and local fields can be explored 
without the extreme dem and on computational complexity involved in multi- 



compartmental modeling (|Linden et al.ll2010l . iLinden et al.ll201ll iProtopapas et al. 



19981 . ISargsyan et al.ll200ll ) by adopting reduced compartment circuits. For 



example, we envisage to extend our recent work which maps fir ing rate model 



(deri ved from LIF models) to population density models (jChizhov et al 



20071 ). but now incorporating our observational DFP model. In addition, 



our framework is analytically amenable and thus can be applied to any lin- 
ear differential equation, for instance, GIF (Gif-sur-Yvette Integrate Fire) 
models, which are improvements to the LIF models and compute more ac- 



2(3 



curately spike activations ( iRudolph-Lilith et al.ll2012l ). Also resonant mem- 
branes (mediated by Ca 2+ and a Ca 2+ -activated K + ionic currents) that 
describe sub-th reshold oscillations and which can be easily expressed by lin- 
ear equations (IMauro et al.l Il970l ) can be incorporated in our derivations. 
We note however that our f ramew ork can be applied to non-linear equations, 
with iHodgkin and Huxleyl ( 119521 ) type activation, but it will fall short from 
explicit and analytical observation equations. 
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